MODULE global_variables
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! Contains all the global variables that will be used in the program !
  !                                                                    !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  !     qphll                                                          !
  !                                                                    !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Number of particles  
!----------------------------------------------------------------------!  
  ! Number of nitrogen molecules
  INTEGER :: n_molecules 
  
  INTEGER,PARAMETER :: n_molecules_max=1000 
  
! Positions
!----------------------------------------------------------------------!
  ! Positions of the centers of mass of the N molec.
  DOUBLE PRECISION,DIMENSION(n_molecules_max) :: rx, ry, rz 
  
! Geometry  
!----------------------------------------------------------------------!
  ! Simulation box length
  DOUBLE PRECISION :: box
! Lennard Jones parameters
  DOUBLE PRECISION :: sig
  DOUBLE PRECISION :: eps
! Physical variables
!----------------------------------------------------------------------!
  ! Temperature (non-dimensional)
  DOUBLE PRECISION :: temp
  ! Energy, etc
  DOUBLE PRECISION :: energy
  ! Block accumulators
  DOUBLE PRECISION :: blkener, blkpress, blkenersq
  ! Virial and pressure
  DOUBLE PRECISION :: pressure
! Parameters
!----------------------------------------------------------------------!
  ! Cutoff distances
  DOUBLE PRECISION :: cutoff_nn, cutoff_nn_square
  ! Very small distances will cause huge energies
  DOUBLE PRECISION, PARAMETER :: huge_number = 10000
  DOUBLE PRECISION :: small_rsq_nn, huge_energy
! Acceptance/Rejections Statistics
!----------------------------------------------------------------------!
  ! Translation statistics
  INTEGER :: acctrans, ntrans
! Monte Carlo parameters
!----------------------------------------------------------------------!
  ! Number of blocks
  INTEGER :: nblocks
  ! Number of steps per block
  INTEGER :: nsteps
  ! Maximum displacement in a MC move
  DOUBLE PRECISION :: maxdisp
  ! Adjust the parameter maxdisp every nadjust blocks
  INTEGER :: nadjust
! Radial distribution function
!----------------------------------------------------------------------!
  INTEGER,PARAMETER :: maxbin = 200
  INTEGER, DIMENSION(maxbin) :: hist
  DOUBLE PRECISION :: delr,constant
  INTEGER :: nsamplegr,totalsamplesgr
! DCD file output
!----------------------------------------------------------------------!
  INTEGER :: dcd_freq
! Others
!----------------------------------------------------------------------!
  ! Dummy variable for random number generator
  INTEGER :: idum
  ! Random number generator
  DOUBLE PRECISION,EXTERNAL :: random
  ! Pi 
  DOUBLE PRECISION,PARAMETER :: Pi = 3.141592d0
  DOUBLE PRECISION,PARAMETER :: Kb = 1.3806d-23  ! J/K
  DOUBLE PRECISION,PARAMETER :: Na = 6.0221d23   ! 1/mol
  DOUBLE PRECISION,PARAMETER :: R = Kb*Na        ! J/(mol K)
  DOUBLE PRECISION,PARAMETER :: Pp = Kb*1.0d+30  ! (J cu.A)/(K cu.M) OR (Pa cu.A)/K
  DOUBLE PRECISION,PARAMETER :: Pb = Pp*1.0d-5   ! (bar cu. A)
END MODULE global_variables
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
PROGRAM nvt
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! Main program                                                       !
  ! Procedures used:                                                   !
  !  - start, runner, finish                                           !
  !                                                                    !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! qphll                                                              !
  !                                                                    !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
  USE global_variables
  IMPLICIT NONE
  CALL start
  CALL runner
  CALL finish
  STOP
END PROGRAM nvt
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE runner
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! Performs 'nblocks' blocks of 'nsteps' Monte Carlo steps            !
  !                                                                    !
  ! Procedures used:                                                   !
  !  - random                                                          !
  !  - translate                                                       !
  !  - total_energy                                                    !
  !  - adjust_mc_parameters                                            !
  !                                                                    !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! qphll                                                              !
  !                                                                    !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
  USE global_variables
  IMPLICIT NONE
  ! Local variables
  INTEGER :: iblock, istep
  DOUBLE PRECISION :: etest, ptest, Utest, Wtest
  ! Calculate the total energy
  CALL total_energy(energy,pressure)
  ! Perform runs over blocks
  DO iblock=1,nblocks
     ! Set block accumulators to zero
     blkener = 0.0d0
     blkenersq = 0.0d0
     blkpress = 0.0d0
     ! Perform runs over steps
     DO istep=1,nsteps
        ! Perform an MC step
        CALL translate
        ! Update block accumulators
        blkener = blkener + energy
        blkpress = blkpress + pressure
        blkenersq = blkenersq + energy**2
        ! End loop over steps
     END DO
        ! Average properties at the end of the block
        blkener = blkener / DBLE(nsteps)/DBLE(n_molecules)
        blkenersq = blkenersq / REAL(nsteps)/DBLE(n_molecules)**2
        blkpress = blkpress / DBLE(nsteps)
    CALL total_energy(energy,pressure)
        ! Write bulk properties to file
        IF(sig == 1.0 .AND. eps == 1.0) THEN
          ! Write results in output file
          WRITE(40,FMT='(I5,F18.6,F24.6,F18.6)') &
               &iblock,blkener,blkenersq,&
               &((n_molecules*temp-blkpress)/box**3)
        ELSE
          ! Write results in output file
          WRITE(40,FMT='(I5,F18.6,F24.6,F18.6)') &
               &iblock,R*blkener,(R**2)*blkenersq,&
               &Pb*((n_molecules*temp-blkpress)/box**3)
        END IF
        ! Adjust MC parameters
        IF (MOD(iblock,nadjust) == 0) CALL adjust_mc_parameters
        ! Sample radial distribution function
        IF (MOD(iblock,nsamplegr) == 0) CALL sample_gr
        ! Write to the DCD file
        CALL write_dcd(iblock)
    
     ! End loop over blocks
  END DO
  ! Close output file
  CLOSE(UNIT=40)
  RETURN
END SUBROUTINE runner
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
REAL FUNCTION random(idum)
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! Generates a random number between 0.0 and 1.0                      !
  ! Modification of ran3, taken from "Numerical Recipes in Fortran"    !
  !                                                                    !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! qphll                                                              !
  !                                                                    !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
  IMPLICIT NONE
  
  ! Local variables
  REAL, PARAMETER :: mbig=4000000., mseed=1618033., mz=0., fac=1./mbig
  INTEGER :: idum,i, ii,  k
  INTEGER, SAVE :: iff, inext, inextp
  REAL :: mj, mk
  REAL, SAVE :: ma(55)
  DATA iff /0/
  ! Initialization
  IF(idum<0.OR.iff==0) THEN
     iff=1
     mj=mseed-IABS(idum)
     mj=MOD(mj, mbig)
     ma(55)=mj
     mk=1
     DO i=1,54
        ii=MOD(21*i,55)
        ma(ii)=mk
        mk=mj-mk
!!modified by lgk on 17 Sept., 2013 
        IF(mk<mz)mk=mk+mbig       
          mj=ma(ii)
!!modified by lgk on 17 Sept., 2013 
      END DO
      DO k=1,4
        DO i=1,55
           ma(i)=ma(i)-ma(1+MOD(i+30,55))
           IF(ma(i)<mz)ma(i)=ma(i)+mbig  
        END DO
     END DO
     inext=0
     inextp=31
     idum=1
  END IF
  ! Start to generate the random number
  inext=inext+1
  IF(inext==56) inext=1
  inextp=inextp+1
  IF(inextp==56) inextp=1
  mj=ma(inext)-ma(inextp)
  IF(mj<mz)mj=mj+mbig 
  ma(inext)=mj
  random=mj*fac
  RETURN
END FUNCTION random
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE  initialize_idum(idum)
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! Initializes the seed of the random number generator                !
  !                                                                    !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! qphll                                                              !
  !                                                                    !
  !                                                                    !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
  IMPLICIT NONE
  ! Local variables
  INTEGER :: idum
  INTEGER, DIMENSION(8) :: val
  idum = -478
!!$  CALL DATE_AND_TIME(VALUES=val)
!!$  idum=-(1000*val(7)+val(8))
  RETURN
END SUBROUTINE initialize_idum

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	SUBROUTINE translate
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! Generates a random displacement of a molecule itest                !
  ! selected randomly                                                  !
  !                                                                    !
  ! Procedures used:                                                   !
  !  - random(idum): random number generator                           !
  !  - energy_calculation                                              !
  !                                                                    !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! qphll                                                              !
  !                                                                    !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
  USE global_variables
  IMPLICIT NONE
  ! Local variables
  INTEGER :: imol
  DOUBLE PRECISION :: rx_old, ry_old, rz_old, Wcr
  DOUBLE PRECISION :: rx_new, ry_new, rz_new
  DOUBLE PRECISION :: vold, vnew,deltav, virialold,virialnew
  LOGICAL :: accepted
  ! Select the particle to be translated
  imol = INT(random(idum) * n_molecules) + 1
  
  ! Read the current position of the molecule itest
  rx_old = rx(imol)
  ry_old = ry(imol)  
  rz_old = rz(imol)
  ! Calculate energy with the old position
  CALL energy_calculation(imol,rx_old,ry_old,rz_old,vold,virialold)
  ! Assign the new position randomly
  ! The maximum displacement is given by maxdisp
  rx_new = rx_old + maxdisp * (2.0d0*random(idum)-1.0d0)
  ry_new = ry_old + maxdisp * (2.0d0*random(idum)-1.0d0)
  rz_new = rz_old + maxdisp * (2.0d0*random(idum)-1.0d0)
  ! Apply periodic boundary conditions
  rx_new = rx_new - DNINT(rx_new/box)*box
  ry_new = ry_new - DNINT(ry_new/box)*box
  rz_new = rz_new - DNINT(rz_new/box)*box
  ! Calculate the energy with the new position
  CALL energy_calculation(imol,rx_new,ry_new,rz_new,vnew,virialnew)
  ! Calculate difference in potential energy
  deltav = vnew - vold
  ! Apply the acceptance criteria
!   ------------------   whl  ----------------------	
	Wcr = random(idum)
  IF (deltav < 0.0d0) THEN
     accepted = .TRUE.
  ELSEIF (deltav >= (huge_energy-1.0d0)) THEN
     accepted = .FALSE.
  ELSEIF ( Wcr < dEXP( -deltav/temp )) THEN
     accepted = .TRUE.
  ELSE
     accepted = .FALSE.
  END IF
!            ----------------------------------------
  ! If the move is accepted...
  IF (accepted) THEN
     
     ! Update the number of accepted translation moves
     acctrans = acctrans + 1
     ! Update positions
     rx(imol) = rx_new
     ry(imol) = ry_new
     rz(imol) = rz_new
     ! Update the value of the energy
     energy = energy + deltav
     ! Update the value of pressure
     pressure = pressure + virialnew - virialold
  END IF
  ! Update the number of translation moves 
  ntrans = ntrans + 1
!WRITE(*,10) imol, vold, vnew, energy, pressure 
!10 FORMAT(I3, 2X F11.6, 2X, F11.6, 2X, F11.6, 2X, F11.6, 2X, F11.6)
  RETURN
END SUBROUTINE translate
SUBROUTINE energy_calculation(itest,rxi,ryi,rzi,potenergy,virial)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! Calculates the potential energy of molecule"i" (test molecule)     !
  ! with all the other molecules/atoms ("j")                           !
  !                                                                    !
  !                                                                    !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! qphll                                                              !
  !                                                                                                                                       !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
  USE global_variables
  IMPLICIT NONE
  ! Variables with argument association
  INTEGER :: itest
  DOUBLE PRECISION :: rxi,ryi,rzi,potenergy,virial
  ! Local variables
  DOUBLE PRECISION :: energy_ij,virial_ij
  DOUBLE PRECISION :: factor
  DOUBLE PRECISION :: rxij,ryij,rzij,rijsq
  INTEGER :: j
  ! Contribution due to interactions with other nitrogen molecules
  !Initialize variables
  potenergy = 0.0d0
  virial = 0.0d0
  ! Start loop over all nitrogen molecules
  DO j=1,n_molecules
     ! avoid self-interactions
     IF (j==itest) CYCLE
     ! calculate vector between centers of mass
     rxij = rx(j) - rxi
     ryij = ry(j) - ryi
     rzij = rz(j) - rzi
     ! apply minimum image convention
     rxij = rxij - DNINT(rxij/box)*box
     ryij = ryij - DNINT(ryij/box)*box
     rzij = rzij - DNINT(rzij/box)*box
     ! Calculate the distance between centers of mass
     rijsq = rxij**2 + ryij**2 + rzij**2
     ! Apply cutoff distance
     IF (rijsq <= cutoff_nn_square) THEN
        ! Calculate the potential due to dispersion interactions
        ! Check if the distance is very small
        IF (rijsq <= small_rsq_nn) THEN
           potenergy = huge_energy
           RETURN
        END IF
        ! Evaluate the L-J potential
        factor = (sig**2/ rijsq)**3
        energy_ij = factor*(factor-1.0d0)
        ! Evaluate the virial 
        virial_ij = factor*(1.0d0-2.0d0*factor)
        ! Add contributions to energy and virial
        potenergy = potenergy + energy_ij
        virial = virial + virial_ij
        ! End the calculations for distance smaller than cutoff
     END IF
     ! End loop over nitrogen molecules
  END DO
  ! Multiply energy and virial by the L-J factor
  potenergy = 4.0d0 * eps * potenergy
  virial = 8.0d0 * eps * virial
END SUBROUTINE energy_calculation
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE total_energy(potenergy,virial)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! Calculates the total energy of the system                          !
  !                                                                    !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! qphll                                                              !
  !                                                                                                                                           !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
  USE global_variables
  IMPLICIT NONE
  ! Variables with argument association
  DOUBLE PRECISION :: potenergy,virial
  ! Local variables
  DOUBLE PRECISION :: rxi,ryi,rzi
  DOUBLE PRECISION :: energy_ij, virial_ij
  DOUBLE PRECISION :: factor
  DOUBLE PRECISION :: rxij,ryij,rzij,rijsq
  INTEGER :: j,imol
  ! Initialize parameters
  potenergy = 0.0d0
  virial = 0.0d0
  ! Start loop over all nitrogen molecules
  DO imol=1,n_molecules-1
     rxi = rx(imol)
     ryi = ry(imol)
     rzi = rz(imol)
     ! Contribution due to interactions with other nitrogen molecules
     ! Start loop over other nitrogen molecules
        DO j=imol+1,n_molecules
           ! calculate vector between centers of mass
           rxij = rx(j) - rxi
           ryij = ry(j) - ryi
           rzij = rz(j) - rzi
           ! apply minimum image convention
           rxij = rxij - DNINT(rxij/box)*box
           ryij = ryij - DNINT(ryij/box)*box
           rzij = rzij - DNINT(rzij/box)*box
           ! Calculate the distance between centers of mass
           rijsq = rxij**2 + ryij**2 + rzij**2
           ! Apply cutoff distance
           IF (rijsq <= cutoff_nn_square) THEN
              ! Calculate the potential due to dispersion interactions
              ! Check if the distance is very small
              IF (rijsq <= small_rsq_nn) THEN
                 energy = huge_energy
                 PRINT*, 'ERROR: Try with a different initial config'
                 STOP
              END IF
              ! Evaluate the L-J potential
              factor = (sig**2/ rijsq)**3
              energy_ij = factor*(factor-1.0d0)
              ! Calculate the virial
              virial_ij = factor*(1.0d0-2.0d0*factor)
              ! Add contributions to energy and pressure
              potenergy = potenergy + energy_ij
              virial = virial + virial_ij
              ! End the calculations for distance smaller than cutoff
           END IF
           ! End loop over nitrogen molecules
        END DO
     ! End loop over all nitrogen molecules
  END DO
  ! Multiply energy and pressure by the L-J factor
  potenergy = 4.0d0* eps * potenergy
  virial = 8.0d0* eps * virial
END SUBROUTINE total_energy
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
SUBROUTINE adjust_mc_parameters
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! Adjusts MC parameter maxdisp                                       !
  !                                                                    !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! qphll                                                              !
  !                                                                                                                                           !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
  USE global_variables
  IMPLICIT NONE
  ! Adjust maxdisp
  IF ( (DBLE(acctrans)/DBLE(ntrans)) > 0.50d0) THEN
     IF (maxdisp < box)  maxdisp = 1.050d0 * maxdisp
  ELSE
     maxdisp = 0.950d0 * maxdisp
  ENDIF
  RETURN
END SUBROUTINE adjust_mc_parameters
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
SUBROUTINE sample_gr
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! Updates the histogram for g(r)                                     !
  !                                                                    !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! qphll                                                              !
  !                                                                                                                                           !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
  USE global_variables
  IMPLICIT NONE
  ! Local variables
  INTEGER :: imol,j
  DOUBLE PRECISION :: rxi,ryi,rzi
  DOUBLE PRECISION :: rxij,ryij,rzij,rijsq,rij
  INTEGER :: bin
  totalsamplesgr = totalsamplesgr + 1
  ! Start loop over all nitrogen molecules
  DO imol=1,n_molecules-1
     rxi = rx(imol)
     ryi = ry(imol)
     rzi = rz(imol)
     ! Contribution due to interactions with other nitrogen molecules
     ! Start loop over other molecules
     DO j=imol+1,n_molecules
        ! calculate vector between centers of mass
        rxij = rx(j) - rxi
        ryij = ry(j) - ryi
        rzij = rz(j) - rzi
        ! apply minimum image convention
        rxij = rxij - DNINT(rxij/box)*box
        ryij = ryij - DNINT(ryij/box)*box
        rzij = rzij - DNINT(rzij/box)*box
        ! Calculate the distance between centers of mass
        rijsq = rxij**2 + ryij**2 + rzij**2
        rij = dSQRT(rijsq)
        IF (rij < box/2) THEN
           bin = INT(rij/delr) + 1
           hist(bin) = hist(bin) + 2
        END IF
        ! End loop over other molecules
     END DO
     ! End loop over molecules
  END DO
  RETURN
END SUBROUTINE sample_gr
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
SUBROUTINE start
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! Sets everything up to start running the MC code                    !
  !                                                                    !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! qphll                                                              !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
  USE global_variables
  IMPLICIT NONE
  ! Local variables
  INTEGER :: imol,bin
  CHARACTER(LEN=5) :: atype
  ! Read the input file
  OPEN(UNIT=10,FILE='inputfile.txt',STATUS='OLD',ACCESS='SEQUENTIAL',&
       &ACTION='READ')
  READ(10,*)
  READ(10,*) temp
  READ(10,*) cutoff_nn
  READ(10,*)
  READ(10,*)
  READ(10,*) sig, eps
  cutoff_nn_square = cutoff_nn * cutoff_nn
  READ(10,*)
  READ(10,*)
  READ(10,*) nblocks
  READ(10,*) nsteps
  READ(10,*) nadjust
  READ(10,*) nsamplegr
  READ(10,*) dcd_freq
  READ(10,*) maxdisp
  
  CLOSE(UNIT=10)
  ! Calculate the very small distances that will cause huge energies
  huge_energy = temp * huge_number
  small_rsq_nn = (0.80d0)**2
  ! Get initial positions of nitrogen molecules
  OPEN(UNIT=20,FILE='initial.xyz',STATUS='OLD',ACCESS='SEQUENTIAL',&
	     &ACTION='READ')
  READ(20,*) n_molecules
  READ(20,*) box

  ! Read positions
  DO imol=1,n_molecules
     READ(20,*) atype,rx(imol),ry(imol),rz(imol)
  END DO
  CLOSE(UNIT=20)
  ! Initialize dummy variable for random number generator
  CALL initialize_idum(idum)
  ! Initialize acc/rej statistics
  acctrans = 0
  ntrans = 0
  ! Initialize radial distribution function
  DO bin=1,maxbin
     hist(bin) = 0
  END DO
  delr = box/(2*maxbin)
  constant=4*pi*DBLE(n_molecules)/(3.0d0*box**3)
  totalsamplesgr = 0

  ! Open the DCD file for writing
  CALL open_dcd
  ! Open file to write results
  OPEN(UNIT=40,FILE='blkaverages',STATUS='UNKNOWN',ACCESS='SEQUENTIAL')
  WRITE(40,*) '#NVT Results (blocks averages)'
  WRITE(40,*) 
  WRITE(40,FMT='(F14.6, '' = Temperature'')') temp
  WRITE(40,FMT='(9X,I5, '' = Number of molecules'')') n_molecules
  WRITE(40,FMT='(F14.6, '' = Volume of the simulation box'')') box**3
  WRITE(40,*) 
  IF(sig == 1.0 .AND. eps == 1.0) THEN
    WRITE(40,*) '#block       U (reduced)       U^2 (reduced)          P (reduced) '
  ELSE
    WRITE(40,*) 'block       U (J/mol)         U^2 (J/mol)^2          P (bar)'
  ENDIF
  
  RETURN
END SUBROUTINE start
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE finish
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! Reports the final results of the MC code                           !
  !                                                                    !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! qphll                                                              !
  !                                                                                                                                           !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
  USE global_variables
  IMPLICIT NONE
  ! Local variables
  INTEGER :: imol,j
  DOUBLE PRECISION :: rlow,rup,nideal
  INTEGER :: bin
  OPEN(UNIT=50,FILE='final.xyz',STATUS='UNKNOWN',ACCESS='SEQUENTIAL')
  WRITE(50,FMT='(1I4,12X)') n_molecules
  WRITE(50,*) box,box,box
  DO imol=1,n_molecules
     WRITE(50,FMT='(A5,3F12.6)')'Ar', rx(imol),ry(imol),rz(imol)
  END DO
  CLOSE(UNIT=50)
  ! Calculate radial distribution function
  OPEN(UNIT=60,FILE='gr',STATUS='UNKNOWN',ACCESS='SEQUENTIAL')
  DO bin=1,maxbin
     rlow = delr*(bin-1)
     rup = rlow + delr
     nideal = constant*(rup**3-rlow**3)
     WRITE(60,FMT='(F8.4,F12.6)') rlow + delr/2 , &
          &hist(bin)/nideal/DBLE(n_molecules)/DBLE(totalsamplesgr)
  END DO
  CLOSE(UNIT=60)
     
  RETURN
END SUBROUTINE finish
  SUBROUTINE open_dcd
    USE global_variables
    IMPLICIT NONE
    ! Internal
    INTEGER :: j,n_frames
    INTEGER :: ntitle
    CHARACTER(LEN=80),DIMENSION(100) :: title
    INTEGER,DIMENSION(20) :: ictrl
    REAL*4 :: delta
    DOUBLE PRECISION,DIMENSION(6) :: pbcs
    ! Initialize delta
    delta = 1.0000
    ! Initialize ictrl
    ictrl = 0
    ! Calculate the number of frames
    n_frames = nblocks/dcd_freq
    ! Set ictrl elements
    ictrl(1) = n_frames
    ictrl(3) = 1
    ictrl(4) = ictrl(1)
    ictrl(11) = 0
    ictrl(20) = 24
    ! Set info for filler junk  
    ntitle = 2
    title(1) = "Created by JDM plugin                     "
    title(2) = "REMARKS Created 10 November, 2008 at 19:24"
    ! Set PBC information
    pbcs(1) = DBLE(box)
    pbcs(2) = 0.0d0
    pbcs(3) = DBLE(box)
    pbcs(4) = 0.0d0
    pbcs(5) = 0.0d0
    pbcs(6) = DBLE(box)
    OPEN(UNIT=5010,FILE='movie.dcd',FORM='UNFORMATTED',STATUS='UNKNOWN')
      ! Write the junk at the begining
      WRITE(5010) "CORD",(ictrl(j),j=1,9),delta, (ictrl(j),j=11,20)
      WRITE(5010) ntitle,(title(j),j=1,ntitle)
      WRITE(5010) n_molecules
  !   WRITE(5010) (pbcs(j),j=1,6)
    RETURN
  END SUBROUTINE open_dcd
  SUBROUTINE write_dcd(iblock)
    USE global_variables
    IMPLICIT NONE
    ! Passed
    INTEGER :: iblock
    ! Internal
    INTEGER :: j, n_frames, iframe
    ! Calculate the number of frames
    n_frames = nblocks/dcd_freq
    ! Calculate the current frame number
    iframe = iblock/dcd_freq
    ! Check if a step should be written at iblock
    IF(MOD(iblock,dcd_freq) == 0 ) THEN
      ! Write the coordinates to file (must be REAL*4)
      WRITE(5010)(REAL(rx(j)),j=1,n_molecules)
      WRITE(5010)(REAL(ry(j)),j=1,n_molecules)
      WRITE(5010)(REAL(rz(j)),j=1,n_molecules)
    END IF
    ! Close the file if it is the last frame to be written
    IF(iframe == n_frames) CLOSE(UNIT=5010)
    RETURN
  END SUBROUTINE write_dcd
